Load packages and configure the grid resolution and source filename:

library(dplyr)
library(ggplot2)
library(dggridR)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(prettygrids)
world <- ne_countries(scale = "medium", returnclass = "sf")

resolution <- 7
filename <- "occurrence_minimal_20200212"

Use the dggridR package to create a discrete global grid, and fetch the grid cell surface area to be included in the map title:

dggs <- dgconstruct(projection = "ISEA", res = resolution, resround = "up")
area <- format(signif(dggetres(dggs)$area_km[resolution + 1], 2), big.mark = ",")

This notebook uses a A CSV file with occurrence records exported from the database database. Alternatively the robis R package can be used, but downloading the entire database can take a while. This downloads the zipped CSV file and aggregates the data into locations with number of records:

if (!file.exists("grouped.temp")) {
  if (!file.exists("data.zip")) {
    download.file(paste0("https://download.obis.org/export/", filename, ".zip"), "data.zip")
  }
  df_raw <- read.csv(unz("data.zip", paste0(filename, ".csv")), stringsAsFactors = FALSE)
  df_grouped <- df_raw %>%
    group_by(decimallongitude = round(decimallongitude, 2), decimallatitude = round(decimallatitude, 2), speciesid) %>%
    summarize(records = n(), species = length(unique(speciesid)))
  save(df_grouped, file = "grouped.temp")
  file.remove("data.zip")
} else {
  load("grouped.temp")
}

Next we assign cell IDs to each location and caluclate statistics per cell:

df_grouped$cell <- dgtransform(dggs, df_grouped$decimallatitude, df_grouped$decimallongitude)
stats <- df_grouped %>%
  group_by(cell) %>%
  summarise(records = sum(records), species = length(unique(speciesid)))

Maps

map_theme <- theme(
  axis.ticks.length = unit(0, "pt"),
  line = element_blank(), rect = element_blank(),
  axis.text = element_blank(), axis.title = element_blank()
)

map_scale_records <- scale_fill_viridis_c(
  option = "plasma",
  trans = "log10",
  name = "Records",
  breaks = c(1, 10, 100, 1000, 10000),
  labels = function(x) format(x, scientific = FALSE),
  limits = c(1, 100000),
  oob = scales::squish
)

map_scale_species <- scale_fill_viridis_c(
  option = "plasma",
  trans = "log10",
  name = "Species",
  labels = function(x) format(x, scientific = FALSE)
)

Pacific centered Robinson projection

Depending on the projection the grids produced by dggridR can cause artefacts in your map. Here we are using the prettygrids::make_grid utility function (still under development) to fix some of these problems.

grid_pac <- make_grid(offset = -180, res = resolution) %>%
  merge(stats, by.x = "cell", by.y = "cell")

crs <- "+proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
bbox <- st_bbox(st_as_sfc(st_bbox(c(xmin = 18, xmax = 248, ymin = -55, ymax = 55), crs = 4326)) %>% st_transform(crs))

Records

ggplot() +
  geom_sf(data = grid_pac %>% st_transform(crs), aes(fill = records), size = 0, color = NA) +
  geom_sf(data = world %>% prepare_grid(offset = -180) %>% st_transform(crs), fill = "#cccccc", color = "#999999", size = 0.1) +
  coord_sf(crs = crs, xlim = c(bbox$xmin, bbox$xmax), ylim = c(bbox$ymin, bbox$ymax)) +
  map_theme +
  map_scale_records +
  labs(
    title = "Number of records",
    subtitle = paste0("Number of records per grid cell of approximately ", area, " km2, Pacific centered Robinson projection. Ocean Biogeographic Information System (OBIS), 2020.")
  )

Species

ggplot() +
  geom_sf(data = grid_pac %>% st_transform(crs), aes(fill = species), size = 0, color = NA) +
  geom_sf(data = world %>% prepare_grid(offset = -180) %>% st_transform(crs), fill = "#cccccc", color = "#999999", size = 0.1) +
  coord_sf(crs = crs, xlim = c(bbox$xmin, bbox$xmax), ylim = c(bbox$ymin, bbox$ymax)) +
  map_theme +
  map_scale_species +
  labs(
    title = "Number of species",
    subtitle = paste0("Number of species per grid cell of approximately ", area, " km2, Pacific centered Robinson projection. Ocean Biogeographic Information System (OBIS), 2020.")
  )

Lambert azimuthal equal-area projection

north <- st_set_crs(st_as_sf(as(raster::extent(-180, 180, 0, 90), "SpatialPolygons")), 4326)
grid_ortho <- make_grid(offset = NA, res = resolution) %>% 
  merge(stats, by.x = "cell", by.y = "cell") %>%
  filter(st_intersects(geometry, north, sparse = FALSE))

crs <- "+proj=laea +lat_0=52 +lon_0=0 +x_0=4321000 +y_0=3210000 +datum=WGS84 +units=m +no_defs"
bbox <- st_sfc(st_polygon(list(matrix(c(-140, 40, -60, 40, 60, 40, 140, 40, -140, 40), ncol = 2, byrow = TRUE))), crs = 4326)
bbox_trans <- bbox %>% st_transform(crs)

Records

ggplot() +
  geom_sf(data = grid_ortho %>% st_transform(crs), aes(fill = records), size = 0, color = NA) +
  geom_sf(data = world %>% st_transform(crs), fill = "#cccccc", color = "#999999", size = 0.1) +
  coord_sf(crs = crs,
           xlim = c(min(st_coordinates(bbox_trans)[,1]), max(st_coordinates(bbox_trans)[,1])),
           ylim = c(min(st_coordinates(bbox_trans)[,2]), max(st_coordinates(bbox_trans)[,2]))) +
  map_theme +
  map_scale_records +
  labs(
    title = "Number of records",
    subtitle = paste0("Number of records per grid cell of approximately ", area, " km2, Lambert azimuthal equal-area projection. Ocean Biogeographic Information System (OBIS), 2020.")
  )

Species

ggplot() +
  geom_sf(data = grid_ortho %>% st_transform(crs), aes(fill = species), size = 0, color = NA) +
  geom_sf(data = world %>% st_transform(crs), fill = "#cccccc", color = "#999999", size = 0.1) +
  coord_sf(crs = crs,
           xlim = c(min(st_coordinates(bbox_trans)[,1]), max(st_coordinates(bbox_trans)[,1])),
           ylim = c(min(st_coordinates(bbox_trans)[,2]), max(st_coordinates(bbox_trans)[,2]))) +
  map_theme +
  map_scale_species +
  labs(
    title = "Number of species",
    subtitle = paste0("Number of species per grid cell of approximately ", area, " km2, Lambert azimuthal equal-area projection. Ocean Biogeographic Information System (OBIS), 2020.")
  )